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Abstract 

Phase space symmetries inherent in the statistical theory of ideal magnetohydrodynamic (MHD) tur- 
bulence are known to be broken dynamically to produce large-scale coherent magnetic structure. Here, 
results of a numerical study of decaying MHD turbulence are presented that show large-scale coherent 
structure also arises and persists in the presence of dissipation. Dynamically broken symmetries in MHD 
turbulence may thus play a fundamental role in the dynamo process. 
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I Introduction 

We examine (both ideal and decaying) incompressible, homogeneous, three-dimensional (3-D) magnetohy- 
drodynamic (MHD) turbulence using Fourier spectral methods. In these methods, physical fields such as 
the turbulent velocity u and magnetic induction b are represented by truncated Fourier series on an evenly 
spaced N 3 computational grid: 

= £ 4 < k > 'W b <x) = ^ T b(k)e“- (1) 

0<|k|<if 0<|k|<if 

The wave vectors k in k-space have integer components between —N/ 2 and +N/2, while the corresponding 
discrete set of position vectors x in x-space have components Xj = 2'Krrij/N (j = 1,2,3; 0 < rrij < N ), 
where N is the number of grid points in each dimension. The ‘isotropic truncation’ radius is K , with 
0 < k < K < N/2 and k = |k|, so that retained coefficients fit within a sphere in k-space; the exact value of 
K is set by de-aliasing requirements 1 . 

Fourier spectral methods place MHD turbulence in a periodic box that can represent either a small 
volume interior to a larger volume of plasma, or a complete mass of plasma within a 3-torus (such as the 
universe). Also, isolated plasma (such as a star) may be thought of as residing in a box large enough so that 
any plasma flows or magnetic fields outside of or on the surface of the box have negligible energy. In this 
later case, a periodic box serves as an approximate surrogate of the more complicated physical situation. 
Thus, Fourier models are of value by providing at least qualitative information about laboratory, geophysical 
and astrophysical turbulence, especially if dissipation is included. 

Statistical studies of such Fourier models were initiated by T. D. Lee, in a prescient paper 2 that demon- 
strated the existence of canonical ensembles based on the ideal invariants in the phase space defined by the 
set of independent Fourier modes. The total energy E is one ideal invariant, and two other invariants were 
soon discovered: Hm, the magnetic helicity 3 and He, the cross helicity 4 . Once computational resources 
were sufficient, basic aspects of the ‘absolute equilibrium ensemble theory’ of 3-D MHD turbulence were 
described 5 . The broken symmetry and non-ergodicity contained in this theory when magnetic and cross 
helicity are non-zero were subsequently recognized 6 ’ 7,8,9 . 

Early computational studies of dissipative, helical MHD turbulence, motivated by the dynamo problem, 
found evidence of an inverse cascade of magnetic helicity from small scales to large ones, with a concomitant 
generation of strong, large-scale magnetic fields 10,11 . The presence of helicity is critical to this process, as 
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studies of non-helical MHD turbulence have found only small-scale dynamo action 12 . Broken symmetry 
and non-ergodicity, in addition to their intrinsic interest, appear relevant to theoretically understanding the 
magnetic dynamo, and this is our focus here. We will build on previous results in ideal MHD turbulence 
and examine their relevance to dissipative MHD turbulence. 

A recent paper 13 examined ideal, homogeneous, incompressible, 3-D MHD turbulence in a rotating frame 
of reference with and without a mean magnetic field. There turn out to be five general cases of incompressible, 
homogeneous 3-D MHD turbulence, categorized by the absence or presence of a mean magnetic field B 0 
and/or an overall system rotation f2 0 ; these cases are shown in Table 1. (Here, the cases are slightly 
reordered.) In Table 1, a new invariant, the ‘parallel helicity’ Hp = He — ctHm is indicated, which appears 
when Q 0 = crB 0 7^ 0 (<r is a non-zero real number). In the recent study, very long-time 32 3 simulations 
again demonstrated that ensemble predictions for ideal MHD turbulence were dynamically broken to produce 
energetic large-scale coherent structures, both with and without overall rotation. 

In the present work, dissipation is introduced through non-zero viscosity v in eq. (2) and magnetic 
diffusivity r\ in eq. (3) (see below). Results from long-time non-rotating and rotating 64 3 simulations are 
presented and large-scale coherent structures, as seen in previous ideal runs, are still present. In effect, this is 
a ‘magnetic dynamo process’ that is inherent to MHD turbulence, a process that requires an initial energetic 
stirring of a magneto- fluid, but does not require any continuous external forcing, such as convection (which 
can, nevertheless, provide a stirring mechanism). A magnetic dynamo may thus be ultimately related to 
dynamically broken symmetry inherent within MHD turbulence. 

II Basic Theory 

The non-dimensional form of the incompressible MHD equations in a rotating frame of reference with constant 
angular velocity f2 0 and a mean (i.e., uniform and constant) magnetic induction B 0 are well known 14 . They 

= V x [u x (a; + 2fl 0 ) + j x (b + B 0 )] + jA7 2 u;, (2) 

= Vx[ux(b + B 0 )]+!)V 2 b. (3) 

Here, the turbulent velocity u and magnetic induction b satisfy V u = V b = 0, and the associated vorticity 
u> and electric current j are lj = V x u and j = V x b). Density does not appear because it equals unity; 
again, v is the kinematic viscosity and 77 is the magnetic diffusivity. 

In terms of the Fourier expansions (1), x-space V —7 ?k in k-space and 

o>(k) = <k x u(k), j(k) = ik x b(k), k • u(k) = k • b(k) = 0. (4) 

There are two independent complex components in a vector Fourier coefficient, for a total of four real and 
imaginary parts for each Fourier mode u(k) and b(k). Since the reality of u(x) requires that u(k) = u*(— k), 
etc., where * denotes complex conjugation, the number of independent modes is one half of the total number 
M ss 47 tA" 3 /3 of non-zero k within the ball k < K. The number of independent k is thus A f = \ Af and the 
dimension of the phase space T of independent real and imaginary Fourier coefficients is Nr = 4A/\ For the 
64 3 grid numerical results to be presented here, Af' = 57656 and Nr = 461248. 

If we set v — 0 in (2) and i) = 0 in (3), we obtain the equations of ideal MHD turbulence. While any 
real flows must have v 7^ 0 and ry 7^ 0, a study of ideal turbulence produces some interesting theoretical 
results that help in understanding aspects of real turbulence, particularly at larger scales of the flow where 
dissipation is minimal. Next, we review ideal MHD turbulence and then discuss real MHD turbulence. 
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Case 

Mean Field 

Rotation 

Invariants 

I 

B o = 0 

Cl 0 = 0 

E, H c , H m 

II 

B o ^0 

Cl 0 =0 

E, H C 

III 

B 0 = 0 

Cl 0 yf 0 

E, H m 

IV 

Bq^O 

Cl 0 — crB 0 

E, H P 

V 

B 0 7 ^ 0 

Cl 0 x B 0 ^ 0 

E 


Table 1: Invariants for ideal MHD turbulence. 


Ill Ideal MHD Turbulence 

First, we define the volume average (i.e., the average value per grid point) of a product <^(x)^>(x) in a 
periodic box B of side length 2ir as 

[H\ = [ 0(x)V’(x)d 3 x = -^3 ( 5 ) 

' n> Jb 0<|k|<if 


Note that the sum over k counts each independent contribution twice, since u(— k) = u*(k). Since all 
functions are periodic, we can use equations (1) through (5), along with integration by parts to derive the 
following relations 13,15 : 


dE 

dt 

= —2 (uCl + r/J) , 

(6) 

dH c 

dt 

= Cl 0 ■ [b x u] - ± {v + rf) [j • u:\ , 

(7) 

dH M 

dt 

II 

CO 

0 

X 
,g , 
1 

(8) 


Above, we have the (volume-averaged) energy E, enstrophy Cl, mean-squared current J, cross helicity He 
and magnetic helicity Hm (here, a is defined by b = V x a with V • a = 0): 

E= i [|u| 2 + |b| 2 ] , fi=i[M 2 ], J=i[|j| 2 ], He = 3 [u • b] , H M = i[ a • b] . (9) 

At this point we note that if n o = erB 0 , i.e., if Cl D and B 0 are non-zero and parallel, then eq. (7) can be 
added to — a times eq. (8) to yield 

= -3 (v + v) [j ’ w] -F 0-77 [j ' b ] , h p =H c -<tH m = § [(u — era) b], (10) 

The recently discovered quantity Hp has been named the ‘parallel helicity’ 13,16 . 

When v = r\ = 0, equations (6), (7), (8) and (10) lead us immediately to the invariant integrals for MHD 
turbulence for various values of B 0 and Cl 0 . The five general cases of incompressible, homogeneous, 3-D 
MHD turbulence appear in Table 1. We see that as a mean field or an overall rotation is imposed, or both, 
the number of invariants drop from three to two to one. 

The probability density function D of the absolute equilibrium ensemble, for Case I in Table 1, has the 
canonical form 5 

D = ^exp(—aE — PHc -jHm), Z = J exp(-aE - f3H c - jH M )dT. ( 11 ) 


In eq. (11), Z is the partition function and T is the phase space defined by the Nr independent Fourier 
coefficients. Using eqs. (5) and (9), the integrand in (11) becomes an explicit quadratic form in the modal 
quantities u(k) and b(k); evaluating the Gaussian integral gives the partition function Z: 


z= n ^ 

0<|k| <K 


(<5 4 — a 2 7 2 /fc 2 ) 2 ’ 


S 2 = a 2 - (3 2 / 4. 


(12) 
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In arriving at the above, it is necessary that a > 0 and 6 2 > a\^\/k > 0. Also, the product is over the 
independent k, i.e., k but not — k. Therefore, we have AT' = \ AT factors Z{ k) in eq. (12). 

Using (11), canonical ensemble predictions for moments of the u(k) and b(k) can be made 5 for Case 
I in Table 1. In what follows, the ensemble prediction of a quantity Q will be denoted by ( Q ) and the 
time-average by Q , where 

(Q) = J Q DdT, 


— i r Tf 

Qdt. 

df Jo 


(13) 


The predicted first-order moments are (u(k)) = ^b(k)y = 0 for all cases. The ensemble predictions for 
second-order moments for Case I in Table 1 are given in the following equations (again, S 2 = a 2 — /? 2 /4): 


(a) (|u(k)| 2 ) 


2 a{5 2 — 7 2 /k 2 ) 
<5 4 — a 2 7 2 /k 2 


(b) (|b(k)| 2 ) 


2 aS 2 

S 4 — a 2 7 2 /k 2 ’ 


(c) (u(k) • b*(k)^ 

(d) (a(k)-b*(k)) 


-(35 2 

5 4 — a 2 "f 2 /k 2 

—2^a 2 /k 2 

5 4 — a 2 7 2 /k 2 ' 


(14) 


In (14), modal expectation values of (a) kinetic energy, (b) magnetic energy, (c) cross helicity, and (d) 
magnetic helicity are given for Case I of Table 1. As for the other cases in Table 1, eq. (14) applies 
as follows: Case II, set 7 = 0; Case III, set (3 = 0; Case IV, set 7 = —a/3, so that eq. (11) becomes 
D — Z~ x exp(— aE — (3Hp) ; and in Case V, set f3 = 7 = 0. 

The modal predictions in (14) can be summed to produce the expectation values ( Ek ), {Em), {He), 
and {Hm) of kinetic energy, magnetic energy, cross helicity, and magnetic helicity, respectively; the total 
energy is, of course, {E) = {Ek) + {Em)- First, since ^|b(k)| 2 ^ > (|u(k)| 2 ) term by term, it is immediately 
apparent that {Em) > {Ek) for all cases in Table 1. Next, using eqs. (5), (9) and (14), the following algebraic 
relations can be derived: 

a{E) + (3 {H c ) + 7 {H m ) = 2 r (r = 

2 a {H C ) +13 {E m ) = 0 (15) 


a({E) - 2 {E m )) - 1 {H m ) = 0. 

(For the 64 3 runs to be described presently, r = 0.43988.) The inhomogeneous set of linear equations (15) 
can easily be solved to yield 


a 


(3 


7 


t {E m ) 

{E m ) {{E) - {E m )) - {H c ) 2 


JHc) 

{E m ) 


a 


2 {E m ) - {E) 
{H m ) a 


(16) 


Equations (14), (15) and (16) will be taken as definitive in that they correct various small errors appearing 
previously 7 > 8 >i5,i3_ 

For the constants of the motion in eqs. (16), we will generally use ( E ) = E, {He) = He and {Hm) = 
Hm, where the time-averages come from numerical simulations. Also, recall that {Em) > {Ek), so that 
2 {Em) — {E)) > 0. Thus, using either eq. (14), (15) or (16), we see that we always have f3Hc < 0 and 
7 Hm < 0, regardless of the signs of He and Hm- These relations will be seen presently to have important 
ramifications. 
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The results presented in (16) indicate that the inverse temperatures a, (3 and 7 are not actually inde- 
pendent, but all depend on one variable parameter: {Em)- Thus, if we wish to find the values of a, (3 and 
7 associated with a particular set of (E), {He) and ( Hm ), we must find the expectation value {Em)- This 
can be done after we define the entropy of the system. 

The entropy S follows from the standard formula S = — {DlnD), where D is given by eq. (11). Using 
eqs. (11), (12) and (15), we have (the summation has A f terms) 

S = 2(r +A/Tn7r) — ^ ln((5 4 — a 2 ^ 2 /k 2 ). (17) 

0<|k|<K 

As discussed above, the proper values of a, (3 and 7 in (17) are determined from eqs. (16) and depend on 
{Em), which is not known beforehand. 

Using results discussed by Khinchin 19 , and later applied to MHD turbulence 6,20 , we define an ‘entropy 
functional’ S'{t) by 

S'{t) = 2(r + A/"ln7r) — ln(<5 /4 — a ,2, y' 2 /k 2 ). (18) 

0<|k|<7C 

In the above equation, S' 2 = o’" — /3' 2 / 4, and a! , (3' , 7' are functions of r in the same way as a, f3 and 7 in 
(16) are functions of {Em)- The functional S'{t) has one global minimum at r = {Em)', this minimum can 
found numerically to yield the entropy S = S' {{Em ))■ 

The results in this section comprise established results which are, nevertheless, not necessarily well known. 
This has motivated the relatively brief and hopefully sufficient presentation given above, as well as the next 
section, where we discuss some novel features due to the pseudoscalar constants of the motion. 

IV Broken Ergodicity 

Now, we consider the effect of the classical symmetry transformations of charge inversion C, spatial reflection 
(or parity ) P and time reversal T. These discrete transformations have played an essential role in quantum 
physics 17 ; as it turns out, they are also essential for understanding homogeneous turbulence 9 . Table 2 shows 
the effect of P, C and T on various quantities in MHD turbulence. It can easily be seen that the equations 
of ideal MHD [(2) and (3) with v = 77 = 0] are invariant under P, C or T . However, the helicities He, 
Hm and Hp change sign under P and/or C , but not T (therefore T does not appear to play a fundamental 
role in ideal turbulence). Also, recall that Hp = He — ctHm, where <7 is defined by fi 0 = <rB 0 ; a is thus a 
pseudoscalar under C but not P or T. 

Since He, Hm and Hp are pseudoscalars under P and/or C, it is not immediately apparent that D 
as given in eq. (11) is invariant under P and C. However, as was mentioned above following eq. (16), the 
helicities and their associated inverse temperatures always satisfy (3 {He) < 0 and 7 {Hm) < 0. Thus, the 
inverse temperatures (3 and 7 behave in the same manner under P, C and T as their associated helicities. 
The fundamental result is that (3 and 7 transform as pseudoscalars under P and/or C , as is explicitly seen 
in eqs. (16). Since (3 and 7 always take on a sign opposite to that of their associated helicities, (3 {He) and 
7 {Hm) are scalars under P or C. 

The novel result 9 of all this is that the ensemble defined by the probability density function D given in 
(11) consists of dynamically unconnected components. When D is used for ensemble averaging, the integral 
over phase space includes these disjoint components that can be identified by ‘set characteristic functions’ 
with the form \C = Hc/\Hc\, Xm = Hm/\Hm\ and \p = Hp/\Hp\ , depending on which case in Table 1 
is under consideration. In other words, the phase space ‘surface of constant energy’ is a union of disjoint 
sets marked by \C = if) Xm = ±1 or xp = ±1- The number of disjoint components for the cases in 
Table 1 is 2 M , where M is the number of invariant helicities listed in each of the five cases (i.e., Case I has 
four components, Cases II, III and IV have two components each, and Case V has one component). The 
Birkhoff-Khinchin theorem 19 , tells us that an ensemble is ergodic if and only if the surface of constant energy 
has no disjoint components. Therefore, in Table 1, only Case V is ergodic, while Cases I, II, III and IV are 
non-ergodic. 

As usual for canonical ensembles, although the expectation value (13) is integrated over all phase space, 
it takes its essential value only from the components of constant P, ±Hc and ±Hm (on these components, 
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I 

P 

C 

T 

u 

+ 

- 

+ 

- 

U) 

+ 

+ 

+ 

- 

b 

+ 

+ 

- 

- 

j, a 

+ 

- 

- 

- 

He 

+ 

- 

- 

+ 

Hm 

+ 

- 

+ 

+ 

H P 

+ 

- 

- 

+ 


Table 2: Effect of P, C, T on the signs of various quantities. 


E = E, He = He and Hm = Hm )• As discussed above, the sign taken by (3 and 7 on each component is 
as required by the conditions (3Hc < 0 and 7 Hm < 0. Thus, under P and C, the scalars (3 He and 7 Hm 
have the same non-positive value on each component, so that the statistical theory is invariant under P and 
C (and T), as are the ideal MHD equations. That the expectation value of the helicities is invariant is also 
shown, for example, by the relation (He) = — din Z/d/3, since both (3 and (He) change sign under P and C . 
This invariance under P and C also requires the expectation values of all odd-orcler moments of the phase 
variables to be zero. Yet, when a numerical simulation is run, the time evolution is only over one of the 
disjoint components, i. e., the theory is non-ergodic. 

We can illustrate this disjointness and non-ergodicity pictorially for Case II of Table 1, where E and He 
are constants of the motion, but not Hm- We first transform the phase space variables (u(k), b(k) ; 0 < 
jkj < K} into the equivalent set |z + (k), z~(k); 0 < |kj < A'}, where z ± (k) = u(k) ±b(k) are the Elsasser 
variables 3 . Since E and He are constant, so are the Elsasser energies E\ and A2: 

E X = E + 2H C = \ [|z+| 2 ], E 2 = E-2H c =i[\z~\ 2 ]. (19) 

Let us assume He > 0, so that E\ > E 2 . The AY-dimensional phase space T is the direct sum of two \ Np- 
dimensional parts T + and T _ , corresponding to the independent transverse parts of the z + (k) and z~(k), 
respectively, with 0 < |k| < K. Since E\ is constant, the z + (k) reside on a hypersphere Si of radius E\ in 
T+, and similarly, the z~(k) reside on a hypersphere S 2 of radius E 2 in T - . The product T12 = Si x S2 is 
represented as the outer torus of Figure 1. Under P or C, z + (k) ^ ±z~(k), while (3 — > —j3 and He — ► — He ■ 
The result is that hypersphere Si of radius E\ is now in T _ , and hypersphere S 2 of radius E 2 is in T + . The 
resulting hypertorus T21 = S2 x Si is represented by the inner torus of Figure 1. If He yf 0, then T\ 2 
and X21 are disjoint and form a nested pair of hypertori, as qualitatively shown in Figure 1. The union 
T = T12 U X21 is the ‘surface of constant energy’ and is essentially disjoint 9 so that the statistical ensemble 
is non-ergodic 19 . (Figure 1 is also appropriate to Case I of Table 1 with the proviso that the phase point 
is further constrained to lie on the intersections of the nested hypertori with the hypersurfaces defined by 

±| H m — •) 

In contrast to an ensemble prediction arrived at by integrating over all disjoint components labeled by 
different signs of He and/or Hm, when a numerical simulation is run (or a physical system evolves), we 
start with either the plus or the minus sign of each helicity, but not both. Thus, the symmetry of the theory 
under P and C is dynamically broken and mean values of the Fourier coefficients with respect to time are 
not necessarily zero, even though this is the ensemble prediction for ideal MHD turbulence. 

Non-ergodicity caused by broken symmetry has been termed ‘broken ergodicity’ 21 . Numerical examples 
of this for the various ideal cases in Table 1 have been presented recently 13 , where it was seen that relatively 
energetic ideal coherent structures arise primarily in Cases I and III (where the invariant magnetic helicity 
Hm had a value of 0.1398). Let us define ‘coherent energy’ as the energy held in the mean values of the 
Fourier modes. At the final time step of these previous runs, the coherent magnetic energy was 12.9% of 
total energy in Cases I and III, while the coherent kinetic energy was 5.66% of total energy in Case I and 
essentially zero for Case III, where He is not invariant. Less energetic ideal coherent structures occurred and 
were also seen in Case IV, where Hp is invariant. Although the statistical ensemble is manifestly non-ergodic 
for Case II, there was no energetic coherent structure in that case because unrestricted motion on either 
hypertori represented in Figure 1 gives each phase variable a time-averaged value of zero. 
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Figure 1: Phase space for ideal MHD turbulence, Case II of Table 1. 


In this paper, we study only Cases I and III numerically, and do so for both ideal and real MHD. The 
reason for this is that only Cases I and III have B 0 = 0, corresponding to which Hm is an ideal invariant. 
The ideal invariance of Hm is evidently required for the appearance of energetic magnetic coherent structure, 
i.e., for strongly broken symmetry. Our object here is to investigate whether or not this ‘magnetic dynamo 
process’ is present in dissipative 3-D MHD turbulence, as it is in ideal cases. The applicability of this work 
is to understanding the large-scale magnetic fields generated by astrophysical and geophysical objects as a 
whole, when these objects are situated in negligible external magnetic fields. 

V Numerical Results 

Here, we present results for real (i.e., decaying) runs, in addition to same-initial condition ideal runs, to see 
if coherent structure still persists in the presence of dissipation. We move from the 32 s grid recently used 13 
to a 64 3 grid, which permits a larger dissipation wavenumber while still allowing for relatively long-time 
simulations on available computer resources. 

A Fourier spectral transform method 1 with 3 ld -order time-integration 22 was used for these numerical 
simulations. Four different 64 3 runs will be discussed. These are the Case I runs: Run la (v = r/ = 0, 
0 o = B 0 = 0) and Run 1 (v = r) = 0.004, = B a = 0), as well as the rotating Case III runs: Run Ira 

(u = r] = 0, fi 0 = 1, B a = 0) and Run lr (y = r\ = 0.004, = 1, B 0 = 0). (Here, ‘a’ stands for ‘absolute 

equilibrium ensemble’, i.e., ideal runs, and ‘r’ stands for ‘rotating’.) Runs 1, la, lr and Ira went from t = 0 
to 200, with time-step At = 0.001. For Runs 1 and lr, the dissipation wavenumber 14 Ko was Ke> « K or 
less, where K 2 = 910. 

A note on time units: The dimensionless simulation times quoted here can be related to ‘eddy turn-over 
time’ Th = 27r/fi 1 / 2 or ‘crossing time’ Te = ‘lit/E 1 / 2 or some similar unit. In the ideal Runs la and Ira, 
equilibrium fi ~ 240 and E ~ 1, yielding Tq, ~ 0.40 and Te — 27r, so that the dimensionless simulation time 
unit falls between Tq and Te- In the dissipative runs, ft ~ 23 and E ~ 1 early in the run and decrease 
to fi ~ E ~ 0.003 at the end; the corresponding Tq and Te thus increase from ~ 1 to ~ 100. These 
characteristic times may be kept in mind, though we will primarily refer only to dimensionless simulation 
time t. 

All these runs had identical initial conditions in which the u(k) and b(k) were initialized so that Ek(}z) ~ 
EM(k) ~ k 4 exp(— 2k 2 /k 2 ), k p = 6, with random phase. At t = 0, in all runs, total kinetic energy Ek and 
magnetic energy Em were equal, with E = Ek + Em = 1, while the helicities were He = 0.34804 and 
Hm = 0.091968. During the runs, the constants of the motion drifted by less than 0.1%. The ‘constant’ 
values for E, He and Hm are taken to be their time-averages over the ideal runs, and (Em) is found by 
minimizing the entropy functional (18), which is then used to find a, (3 and 7 through eq. (16). The results 
are given in Table 3. 
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Run 

E 

H C 

Hm 

{Em) 

a 

/? 

7 

la 

1.0003 

0.34794 

0.091969 

0.52605 

1.80134 

-2.38283 

-1.01321 

Ira 

1.0006 

0.00000 

0.091968 

0.54627 

0.96825 

0.00000 

-0.96810 


Table 3: Parameters for Runs la and Ira. 


The appearance of coherent structure can be examined directly using phase diagrams , i.e., plots of the real 
versus the imaginary parts of a Fourier coefficient, which are essentially the projections of the ./Vp-dimensional 
phase trajectory onto 2-D planes. Phase diagrams for w y (k), w z (k), b y (k) and 6 z (k) are presented in Figures 
2 and 3 for k = (1, 0, 0), and in Figures 4 and 5 for k = (2, 0, 0); due to eq. (4), <D x (k) = 6 x (k) = 0 whenever 
k = (fc, 0, 0). (Similar phase diagrams were seen for other k.) Equating ‘non-zero mean value’ with ‘coherent 
structure’, Figure 2 shows that there are both kinetic and magnetic structures for Case I, while Figure 3 
shows that there is only magnetic structure for Case III. Additionally, as Figure 6 shows, it is evident that 
coherent structure exists and is manifested primarily in the k = 1 modes (as discussed below). 

The average energy per x-space grid point is equal to one during the ideal runs, and decreases from 
an initial value of one for the dissipative runs. The Fourier transforms (1) are defined so that the average 
corresponding energy per k, with 0 < |kj < K, in k-space is r^ 1 = 2.2733, where r is defined in (15). The 
expectation values (14) predict that the k = 1 modes will have the largest energies, since 5 A — a 2 y 2 //c 2 is 

smallest for k = 1: for Run la, (lu(k)l 2 ) 1 ^ 2 = 58.8 and ^ |b(k)| 2 ^ = 89.0, while for Run Ira, (lu(k)l 2 ) 1 ^ 2 = 

1.43 and ^|b(k)| 2 ^ = 88.9. When k = 2 (fc 2 = 4), we have <( | u(k) | 2 ^ 1//2 = 1.50 and ^jb(k)| 2 ^ = 1.62, 

while for Run Ira, (|u(k)| 2 ) " = 1.43 and (|b(k)| 2 ) = 1.66. There are four independent parts for 

each of the (|u(k)| 2 ) and <^jb(k)| 2 ^, and while the ensemble prediction is for average equality amongst the 
respective parts, for k = 1 this equality is broken because the modes have large mean values, as opposed 
to the ensemble prediction of zero first-order moments. The values of the modal magnitudes given here are 
reflected in the ranges of the various axes in Figures 2, 3, 4 and 5. 

Ensemble predictions that Fourier modes such as b(k) are zero-mean random variables does not match 
the time behavior, particularly for k = 1, and this is shown in Figures 2 to 5. Non-ergodicity (due to the 
broken symmetry) is indicated by the non-zero mean values that are very large compared to the associated 
standard deviations; in the very-long time 32 3 ideal runs recently reported 13 , this ratio was as high as 18, 
and is similar here. After a brief initial transient period, the phase point arrives at an attractor: for Runs 1 
and la, Figure 2 shows that the ideal trajectory arrives at its attractor at t ss 70, while the real trajectory 
peaks at t « 52 (the phase path ends at t = 200). 

Although the location of the ideal attractor depends on initial conditions, the equations of motion and 
statistical theory are invariant under PCT, as well as under galilean transformations (x — » x + a) and 
rotation of the coordinate system. Thus, the attractor can be placed in various equivalent locations and, in 
effect, an ensemble prediction does just this, leading to a prediction of zero mean values for the first-order 
moments. Also, please remember that non-ergodicity results from disjointness 19 , and exists whether or not 
non-zero mean values of first-order moments occur. 

It is evident that the ideal and real modes in Figures 2-5 track each other fairly closely initially, and are 
qualitatively similar over the whole trajectory. This indicates the pertinence of ideal results to real turbulent 
flows, particularly at low values of k. Here, this is evident on 64 3 grids, where v = r\ = 0.004, and should be 
far more evident when larger grids are used, with correspondingly smaller v and g. However, there is always 
a trade-off between grid-size and run-time, so that small grid-sizes are needed if long run-times are desired, 
for example, in statistical studies. 

Also note that Figures 2 to 5 show what appears to be a general result, termed ‘depression of nonlinear- 
ity’ 23 , i.e., the tendency to a force- free state 24 . In eqs. (2) and (3), the nonlinear effects are due to u x u, 
j x b and u x b; these vector cross products become small when u ~ ay j ~ b and u ~ b. For modes 
with k = (fc, 0,0), this translates to w y (k) ~ ±iw z (k), b v ( k) ~ ±i& z (k) and w y (k) ~ ±i6 z (k), respectively. 
In other words, depression of nonlinearity manifests itself for k = (fc,0,0) modes as a ±7r/2 phase shift of 
ujy(h) and h y (k) with regard to u; z (k) and 6 z (k), respectively. If u ~ b, we also expect (k) ~ ±6 y (k) and 
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uj z (k) ~ ±6-(k). These relationships are clearly seen in Figures 2-5. (Since we are dealing with a statistical 
ensemble, we may also view ‘depression of nonlinearity’ as an application of Le Chatelier’s principle 25 .) As 
nonlinearity is depressed, characteristic times associated with nonlinear evolution are increased, essentially 
reducing the usefulness of qualitative temporal measures such as the previously discussed Tq and Te- 

In Figure 6, the time evolution of energy in modes with k 2 = 1, 2, 3 and 4, as well as the total energy for 
all modes, for the various Runs are shown. In particular, the transfer of energy between shells of different 
k is apparent; initially the spectrum was peaked at k 2 ss 6, and Figure 6 shows that it follows an ‘inverse 
cascade’ 5 to smaller and smaller values of k. In the ideal Runs la and Ira, Figure 6 (a) and (c) show that 
equilibrium is almost established by t = 200. In the decaying Runs 1 and lr, Figure 6 (b) and (d) also show 
an inverse cascade, even though the total energy is monotonically decaying: For Run 1, the energy in the 
k 2 = 4, 3, 2 and 1 modes peaks at t « 2, 10, 20 and 50, respectively, while for Run lr, the behavior is similar, 
except that the k 2 = 1 energy profile is flatter and peaks slightly at t ~ 70. Additionally, the k 2 = 2 modal 
energy appears to be dominant for a longer time in Run Ira and lr, as Figure 6 (b) and (d) show. These 
rotating runs also dissipate energy more slowly than the non-rotating runs; this decrease in dissipation for 
rotating turbulence has been previously noted 26 . Overall, Figure 6 shows that the inverse cascade, coupled 
with the attractor evident in Figures 2-5, indicates the growth of a coherent structure or magnetic dynamo. 

Figure 6 indicates that ideal Runs la and Ira have much of their final energy held in the k = 1 modes. 
In fact, most of this k = 1 energy is in the associated coherent structures. In Run la, the coherent magnetic 
energy is 7.83% and coherent kinetic energy is 3.59% of total energy at t = 200, while in Run Ira, the 
coherent magnetic energy is 6.26% and coherent kinetic energy ~ 0% of total energy at t = 200. These 
percentages are similar to what was seen previously 13 for Case I and III runs. In regard to the decaying 
Runs 1 and lr, Figures 2 and 3 show that the k = 1 modes behave in a very coherent manner, while Figure 
6 shows that the energy in dissipative Runs 1 and lr is essentially all in the k = 1 modes at t = 200. Thus, 
the dissipative runs appear to evolve so that almost 100% of their energy in a coherent, albeit decaying, 
structure. 

Please note that Figures 6 (a) and (c) show that the ideal turbulence simulations are not quite yet in 
equilibrium. Using eqs. (14), the energy in the k 2 ■= 1 modes is predicted to be 0.1302 for Run la and 0.09051 
for Run Ira. In comparison, the (fluctuating) data associated with Figures 6 (a) and (c) have reasonably close 
values of 0.1296 and 0.08774, respectively, at t = 200. However, for k 2 = 2, the corresponding equilibrium 
values for Runs la and Ira are 0.0001553 and 0.0001418, while the values at t = 200 are 0.001107 and 
0.003071, respectively. Thus, while the k = 1 modes have essentially arrived at their equilibrium values, for 
k > 1 this requires some more time, as is indicated in the general downward slope of the k 2 = 2, 3 and 
4 modes in Figures 6 (a) and (c). However, the attainment of complete equilibrium is not our goal here, 
as it was in recent 32 3 simulations 13 , where the run-time went from t = 0 to 2000, so that an equilibrium 
spectrum for all k was achieved. 

Next, we consider the evolution of the cross helicity He and the magnetic helicity Hm for the different 
runs. In Figure 7, corresponding to the non-rotating Runs 1 and la, we see that total (i.e., summed over 
all modes) He and Hm are conserved for Run la, as they should be. For the ideal Run la, the difference 
in He and Hm is that the contribution of the k = 1 modes is only about 15% to the total He, while 
the contribution of the k = 1 modes to Hm is essentially 100%. This is the inverse cascade of magnetic 
helicity observed long ago 5 . For the decaying Run 1, Figure 7 shows that all of the value of He and Hm is 
concentrated in the k = 1 modes at the end of the runs. Furthermore, we note that the cusps in Figure 7 
are related to a change in sign of the k = 1 parts of He and Hm, in that these are negative before the cusp 
and positive afterwards. 

In the rotating Runs lr and Ira, He is not an ideal invariant, while Hm remains so. This is clearly seen in 
Figure 8 (a), where the multitude of cusps indicate that the total He and its k =1 part are both fluctuating 
about zero. Figure 8 (b), however, shows that Hm in the rotating Runs lr and la behave very similarly to 
Hm in the non-rotating runs, by comparison with Figure 7 (b). In all runs, ideal and real, magnetic helicity 
concentrates itself in the largest scale ( k = 1) available. Thus, the ideal invariance of Hm, which is the 
common feature between Cases I and III of Table 1, appears essential for the creation of coherent structure 
in both ideal and dissipative MHD turbulence, with or without rotation. 

In fact, the concentration of magnetic helicity is more important in the decaying runs, than in the ideal 
ones, because the ‘normalized helicities’ actually grow in value. Here, the normalization of cross and magnetic 
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Figure 9: Normalized helicities for (a) Runs 1 and la; (b) Runs lr and Ira. 


helicities will be defined in the following manner: 

h c = h M = Hm 1/2 • (20) 

(- E K E M ) 1/2 {E m A) 1 ' 2 

Above, we use Ek = \ [w 2 ] , Em = § [ b 2 ] and A = \ [a 2 ] . 

In Figure 9 (a), the normalized helicities for the ideal Runs la and Ira are presented, while in Figure 9 
(b), the normalized helicities for the dissipative Runs 1 and lr are given. Although dissipation leads to overall 
energy decay, it is clear in Figure 9 that if He and/or Hm are invariant in the ideal runs, they decay less 
quickly than energy in the corresponding dissipative runs. This is the process of ‘selective decay’ 27,28,29,30 , 
by which states of dynamically aligned or anti-aligned u, ui and b arise. 

In fact, ‘dynamical alignment’ and ‘selective decay’ seem equivalent to ‘depression of nonlinearity’. Broken 
symmetry, on the other hand, allows for large mean values of the k = 1 modes. Coupled with ‘selective decay’, 
it leads, as we have seen, to relatively energetic coherent structures. In particular, it leads to large coherent 
magnetic structures in Cases I and III of Table 1. 

VI Conclusion 

Although the simulations discussed were on a relatively small grid of 64 3 points, this allowed for runs from 
t = 0 to 200 (there is always this trade-off in numerical simulations). In the ideal Run la, broken symmetry 
and non-ergodicity (‘broken ergodicity’ 21 ) leads to coherent kinetic and magnetic structures similar to that 
found in recent 32 3 runs 13 that evolved from t = 0 to 2000. In the present work, the real Run 1 had a k = 1 
modal trajectory close to that of ideal Run la, until about t « 50. Similar behavior was observed in the 
rotating Runs Ira and lr, with the difference that only magnetic coherent structure developed. 

These structures manifest themselves at the lowest wave numbers k = 1, that is, at the largest length 
scales of the system, where eqs. (14) predict the largest ideal modal energies. These structures are due 
to ideal invariants, specifically to the pseudoscalars He and Hm , which enter into eqs. (14) through the 
pseudoscalar inverse temperatures (3 and 7 via eqs. (16). As we have seen, the ensemble phase space 
has disjoint components identified by plus and minus signs of the invariant helicities, which allows for the 
symmetry of the phase space to be maintained under the discrete transformations C and P. 

Although ensemble averages are taken over all components in phase space, dynamical evolution is confined 
to only one component, thereby breaking the symmetry inherent in the phase space and ensuring the canonical 
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ensemble is not ergodic. This non-ergodicity is often readily apparent in the non-equivalence of ensemble 
and time averages; however, there are cases where the dynamical system is non-ergodic but where ensemble 
and time averages happen to agree, at least for first and second-order moments. Perhaps this non-ergodicity 
can be seen in the study of third or higher-order moments, or through some other measure. 

We expect that lower dissipation runs on larger grids will produce even more slowly decaying, longer- 
lived and energetic magnetic dynamos. These larger grids will also allow for a clearer view of modal energy 
transfer and modal contribution to emerging coherent structures. Of course, grid sizes cannot be too large if 
very long simulations times are desired. Studies of forced turbulence should also prove interesting, but care 
must be taken so as not to introduce artifacts due only to the particular forcing method chosen. (Boussinesq 
convection is a natural method of forcing, but this moves us away from the five cases listed in Table 1, since 
a new physical field, temperature variation, is introduced.) 

In addition to continuing the study of homogeneous turbulence using Fourier series, investigations with 
other spectral methods, for example using spherical harmonic expansions, may also be very informative. 
It is a straightforward matter to keep track of the means, variances, and higher-order moments of the 
expansion coefficients, whether the spectral method of choice uses sines and cosines, or spherical harmonics. 
Spherical harmonic methods ( e.g ., with Chebyshev expansions in the radial direction) are not uncommon, 
but statistical studies of the associated canonical ensemble of coefficients appear to be absent. 

In summary, what we have seen here using Fourier spectral methods and relatively small grids is the 
persistence in dissipative flows of the effects of broken symmetry and non-ergodicity that were first seen 
in ideal simulations. The result of this ‘broken ergodicity’ in MHD turbulence is that coherent structure 
appears to grow out of initially random and unstructured initial conditions. The most energetic structures 
occur in Cases I and III of Table 1, and these seem to be primarily magnetic in nature. This ‘magnetic 
dynamo process’ is due to the topology of phase space, which is a collection of disjoint sets. Rotation does 
not itself cause this, and the strength of the magnetic structures does not appear to be greatly affected 
by rotation (it is only the kinetic part of the structure that disappears when rotation is present). In fact, 
studies of magnetic stars indicate only a ‘loose . . . correlation between magnetic flux and angular rotational 
velocity’ 31 . 

Of course, the solar interior and the earth’s outer core are not incompressible, homogeneous magneto- 
fluids, so that compressibility and other effects must eventually be taken into account. Here, for now, we must 
content ourselves with some interesting and suggestive results found in the investigation of incompressible, 
homogeneous MHD turbulence. These results probably follow more from the presence of pseudoscalar helical 
invariants in MHD turbulence, and less from the theoretical approximations employed or numerical methods 
used. This essential feature may help these results retain their pertinence when more detailed physical 
models and numerical methods are implemented. 
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